home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_3.6 / subsampleNI / subsampleNI.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  13.1 KB  |  374 lines
  1. /* 
  2.  * subsampleNI.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /* subsampleNI: program performs image size reduction, in particular
  10.  * enabling non-integer subsampling
  11.  */
  12.  
  13. #define DFLTPCT 50.0            /* default subsample percentage */
  14.  
  15. #define G_THRESHPCT_M 100.0     /* default gray thresh pct for matched fltr */
  16. #define G_THRESHPCT_R 100.0     /* default gray thresh pct for rect fltr */
  17. #define G_THRESHPCT_G 60.0      /* default gray thresh pct for Gaussian fltr */
  18.  
  19. #define MAXCOEF 100.0           /* maximum filter coefficient value */
  20. #define RESOLUTION 8            /* subpixel res. of filter */
  21. #define FLTR_VAR_NUM 2.0        /* num. of var.s at max width of fltr */
  22.  
  23. #define MAX_BITS_QUANT    8        /* max. bits of gray quantization */
  24.  
  25. #define ABSF(A) ((A > 0.0) ? A:-(A))  /* floating absolute value */
  26.  
  27. #include <stdio.h>
  28. #include <stdlib.h>
  29. #include <string.h>
  30. #include <images.h>
  31. #include <tiffimage.h>
  32. #include <math.h>
  33. extern void print_sos_lic ();
  34.  
  35. int shrinkGG (unsigned char **, long, long, unsigned char **, long *, long *,
  36.               double, long, double, double);
  37. long usage (short);
  38. long input (int, char **, double *, long *, long *, long *, double *,
  39.             double *, short *, short *);
  40.  
  41. int argc;
  42. char *argv[];
  43. double *pct;                    /* subsample percentage */
  44. long *wSize, *hSize;            /* subsampled image size: width, height */
  45. long *fltrSize;                 /* filter sidelength */
  46. double *fltrVarNum;             /* num. of var.s at max width of fltr */
  47. double *threshG;                /* gray-scale threshold scale factor */
  48. short *qfFlag;                  /* quick filter if =1, or default if 0 */
  49. short *gfFlag;                  /* Gaussian filter if =1, or default if 0 */
  50.  
  51. main (
  52.        int argc,
  53.        char *argv[])
  54. {
  55.   Image *imgI, *imgO;           /* I/O image structures */
  56.   unsigned char **imgIn, **imgGray;  /* input and gray images */
  57.   long widthIn, heightIn;       /* image size */
  58.   long widthOut, heightOut;
  59.   long fltrSize;                /* sidelength of filter */
  60.   long xSize, ySize;            /* subsampled image size */
  61.   double pct,                   /* percentage of row/cols to remain */
  62.     rate,                       /* subsample rate */
  63.     xRate, yRate,               /* rate to reach desired size */
  64.     threshG,                    /* gray-scale threshold scale factor */
  65.     fltrVarNum,                 /* num of var.s at max width of fltr */
  66.     fltrVar;                    /* variance of Gaussian filter [pix] */
  67.   short qfFlag;                 /* quick filter if =1, or dflt */
  68.   short gfFlag;                 /* quick filter if =1, or dflt */
  69.  
  70.   if (input (argc, argv, &pct, &xSize, &ySize, &fltrSize,
  71.              &fltrVarNum, &threshG, &qfFlag, &gfFlag) < 0)
  72.     return (-1);
  73.  
  74. /* allocate input image memory */
  75.   imgI = ImageIn (argv[1]);
  76.   imgIn = ImageGetPtr (imgI);
  77.   heightIn = ImageGetHeight (imgI);
  78.   widthIn = ImageGetWidth (imgI);
  79.  
  80. /* set threshold contrast factors depending on slow/quick */
  81.   if (threshG == -1.0) {
  82.     if (qfFlag == 1)
  83.       threshG = G_THRESHPCT_R;
  84.     else if (gfFlag == 1)
  85.       threshG = G_THRESHPCT_G;
  86.     else
  87.       threshG = G_THRESHPCT_M;
  88.   }
  89.   threshG = threshG / 100.0;
  90.  
  91. /* determine subsampling rate, filter parameters */
  92.   if (xSize == -1 && ySize == -1)
  93.     rate = 100.0 / pct;
  94.   else {
  95.     xRate = (double) widthIn / (double) xSize;
  96.     yRate = (double) heightIn / (double) ySize;
  97.     if (xRate > yRate)
  98.       rate = xRate;
  99.     else
  100.       rate = yRate;
  101.     if (rate < 1.0)
  102.       rate = 1.0;
  103.     pct = 100.0 / rate;
  104.   }
  105.   if (fltrSize == 0)
  106.     fltrSize = 2 * (long) ((100.0 / pct - 1.0) / 2.0 + 0.9999) + 1;
  107.   if (fltrSize % 2 == 0)
  108.     fltrSize += 1;
  109.   fltrVar = ((double) (fltrSize - 1)) / (2.0 * fltrVarNum);
  110.   if (fltrVar == 0.0)
  111.     fltrVar = 1.0;
  112.  
  113. /* determine output size and allocate output image */
  114.   heightOut = (long) ((heightIn - fltrSize) / rate) + 1;
  115.   widthOut = (long) ((widthIn - fltrSize) / rate) + 1;
  116.   imgO = ImageAlloc (heightOut, widthOut, 8);
  117.   imgGray = ImageGetPtr (imgO);
  118.   printf ("Image input size = %dx%d, filter size = %dx%d\n",
  119.           widthIn, heightIn, fltrSize, fltrSize);
  120.  
  121. /* perform subsampling */
  122.   shrinkGG (imgIn, widthIn, heightIn, imgGray, &widthOut, &heightOut,
  123.             pct, fltrSize, threshG, fltrVar);
  124.  
  125. /* write output image to file */
  126.   ImageOut (argv[2], imgO);
  127.  
  128.   printf ("Image output size = %dx%d\n", widthOut, heightOut);
  129.  
  130.   return (0);
  131. }
  132.  
  133. /* SHRINKGG: shrink program for gray-scale input image uses Gaussian filter
  134.  */
  135.  
  136. int
  137. shrinkGG (imgIn, widthIn, heightIn, imgOut, widthOut, heightOut,
  138.           pct, fltrSize, threshG, fltrVar)
  139.      unsigned char **imgIn, **imgOut;  /* input and output images */
  140.      long widthIn, heightIn;
  141.      long *widthOut, *heightOut;
  142.      double pct;                /* percentage of row/cols to remain */
  143.      long fltrSize;             /* sidelength of filter */
  144.      double threshG;            /* gray-scale threshold scale factor */
  145.      double fltrVar;            /* variance of Gaussian filter [pix] */
  146. {
  147.   long halfFltr,                /* half filter sidelength */
  148.     xOut, yOut;                 /* output image coord.s */
  149.  
  150.   long temp;                    /* value of pixel before normaliz. */
  151.   long xMinW, xMaxW,            /* min and max window coord.s */
  152.     yMaxW, yMinW;
  153.   double xIn, yIn;              /* input image subpixel coord.s */
  154.   long xInMid, yInMid;          /* integer middle of fltr on input */
  155.   double rate;                  /* subsample rate */
  156.   long **coef;                  /* fltr coefs for (res * size) dists */
  157.   long halfCoef,                /* no. coef.s in half filter */
  158.     halfCoefPHalf;              /* no. coef.s in half plus half fltr */
  159.   double exponent;              /* exponent of Gaussian for filter */
  160.   long sum,                     /* result of filter conv. at pt */
  161.     sumCoef;                    /* sum of coef.s applied at point */
  162.   struct point end;             /* bottom right of image */
  163.   double endYMHalf, endXMHalf;  /* last x,y coef, minus half filter */
  164.   double distance;              /* dist * res. on Gaussian filter */
  165.   long distX, distY;            /* x,y dist.s within filter mask */
  166.   double xInR, yInR;            /* yIn,xIn scaled by RESOLUTION */
  167.   long i, j;                    /* increment */
  168.   long iR, jR;                  /* increment scaled by RESOLUTION */
  169.   long xMinWR, yMinWR;          /* min. of mask on scaled coord.s */
  170.   double rateR;                 /* scaled rate */
  171.   double sqrt (), exp (), log2 ();
  172.   long nbrSum;                  /* sum of coef.s for ONs in nbrhd */
  173.  
  174. /* compute filter coef.s */
  175.   halfFltr = (fltrSize - 1) / 2;
  176.   halfCoef = halfFltr * RESOLUTION;
  177.   halfCoefPHalf = halfFltr * RESOLUTION + (RESOLUTION + 1) / 2;
  178.  
  179.   if ((coef = (long **) calloc (halfCoefPHalf + 1, sizeof (long *))) == NULL)
  180.       exit (1);
  181.   for (i = 0; i <= halfCoefPHalf; i++)
  182.     if ((coef[i] = (long *) calloc (halfCoefPHalf + 1, sizeof (long))) == NULL)
  183.         exit (1);
  184.  
  185.   for (i = 0; i <= halfCoefPHalf; i++) {
  186.     exponent = ((double) (i)) / (RESOLUTION * fltrVar);
  187.     coef[0][i] = coef[i][0]
  188.       = (long) (MAXCOEF * exp (-exponent * exponent / 2.0) + 0.5);
  189.   }
  190.   for (j = 1; j <= halfCoefPHalf; j++) {
  191.     for (i = 1; i <= halfCoefPHalf; i++) {
  192.       distance = sqrt ((double) (i * i + j * j));
  193.       exponent = distance / (RESOLUTION * fltrVar);
  194.       coef[j][i] = (long) (MAXCOEF * exp (-exponent * exponent / 2.0) + 0.5);
  195.     }
  196.   }
  197.  
  198. /* filter and subsample */
  199.   end.x = widthIn - 1 - halfFltr;
  200.   end.y = heightIn - 1 - halfFltr;
  201.   endYMHalf = (double) (end.y - halfFltr);
  202.   endXMHalf = (double) (end.x - halfFltr);
  203.  
  204.   yInR = halfFltr * RESOLUTION;
  205.   rate = 100.0 / pct;
  206.   rateR = rate * RESOLUTION;
  207.  
  208.   for (yIn = halfFltr, yOut = 0; yIn < endYMHalf; yIn += rate, yInR += rateR) {
  209.     yInMid = (long) (yIn + 0.5);
  210.     yMinW = yInMid - halfFltr;
  211.     yMaxW = yInMid + halfFltr;
  212.     yMinWR = yMinW * RESOLUTION;
  213.     xInR = halfFltr * RESOLUTION;
  214.     for (xIn = halfFltr, xOut = 0; xIn < endXMHalf; xIn += rate, xInR += rateR) {
  215.       xInMid = (long) (xIn + 0.5);
  216.       xMinW = xInMid - halfFltr;
  217.       xMaxW = xInMid + halfFltr;
  218.       xMinWR = xMinW * RESOLUTION;
  219.       sum = sumCoef = nbrSum = 0;
  220.       for (j = yMinW, jR = yMinWR; j <= yMaxW; j++, jR += RESOLUTION) {
  221.         distY = (long) (ABSF (yInR - jR) + 0.5);
  222.         for (i = xMinW, iR = xMinWR; i <= xMaxW; i++, iR += RESOLUTION) {
  223.           distX = (long) (ABSF (xInR - iR) + 0.5);
  224.           if (imgIn[j][i] != 0)
  225.             sum += coef[distY][distX] * imgIn[j][i];
  226.           sumCoef += coef[distY][distX];
  227.         }
  228.       }
  229.       if (sum != 0) {
  230.         temp = (long) ((double) sum / (threshG * sumCoef) + 0.5);
  231.         if (temp > 255)
  232.           temp = 255;
  233.         imgOut[yOut][xOut++] = (unsigned char) temp;
  234.       }
  235.       else
  236.         imgOut[yOut][xOut++] = 0;
  237.     }
  238.     yOut++;
  239.   }
  240.   *heightOut = yOut;
  241.  
  242.   return (0);
  243. }
  244.  
  245.  
  246. /* USAGE: function gives instructions on usage of program
  247.  *   usage: usage (flag)
  248.  *  When flag is 1, the long message is given, 0 gives short.
  249.  */
  250.  
  251. long
  252. usage (flag)
  253.      short flag;                /* flag =1 for long message; =0 for short message */
  254. {
  255.  
  256. /* print short usage message or long */
  257.   printf ("usage: subsampleNI in.img out.img\n");
  258.   printf ("\t\t[-p PCT <%4.1f%%>] OR ([-h HEIGHT] AND/OR [-w WIDTH])\n",
  259.           DFLTPCT);
  260.   printf ("\t\t[-tg GRAY CONTRAST <%3.0f%%>]\n", G_THRESHPCT_M);
  261.   printf ("\t\t[-f FLTRSIZE]\n");
  262.   printf ("\t\t[-qf QUICK FILTER (rectangular filter)]\n");
  263.   printf ("\t\t[-gf GAUSSIAN FILTER (instead of default matched filter]\n");
  264.   printf ("\t\t[-v VARIANCE MULTIPLE (for Gaussian) <%2.0f>]\n", FLTR_VAR_NUM);
  265.   if (flag == 0)
  266.     return (-1);
  267.  
  268.   printf ("\nSUBSAMPLENI filters and subsamples to chosen size.\n");
  269.   printf ("\nFLAGS:\n");
  270.   printf ("\t-p PCT  percentage size of image sidelengths to retain.\n");
  271.   printf ("\t-h HEIGHT -w WIDTH size of image sidelengths to retain.\n");
  272.   printf ("\t\tPCT can be used alone OR (HEIGHT AND/OR WIDTH). When both\n");
  273.   printf ("\t\tHEIGHT and WIDTH are used, the smaller reduction is made.\n");
  274.   printf ("\t-tg GRAY CONTRAST PCT smaller number gives darker image.\n");
  275.   printf ("\t\tSetting -tg to 100%%, turns off constrasting.\n");
  276.   printf ("\t-f FILTER SIZE is square mask size low-pass filter (>= 3).\n");
  277.   printf ("\t-qf QUICK FILTER - quicker rectangular filter instead of\n");
  278.   printf ("\t\tdefault matched filter.\n");
  279.   printf ("\t\tIf filter size is set to 1, this is the quickest and\n");
  280.   printf ("\t\tugliest result.\n");
  281.   printf ("\t\tNote that this filter is used automatically instead\n");
  282.   printf ("\t\tof default for filter sizes larger than 5x5.\n");
  283.   printf ("\t-gf GAUSSIAN FILTER - symmetric Gaussian filter instead of\n");
  284.   printf ("\t\tdefault matched filter.\n");
  285.   printf ("\t-v VARIANCE MULTIPLE is no. of variances of Gaussian filter\n");
  286.   printf ("\t\tat maximum radius in filter window (0 gives uniform wting).\n");
  287.   printf ("\t-L: print Software License for this module\n");
  288.  
  289.   return (-1);
  290. }
  291.  
  292.  
  293. /* INPUT: function reads input parameters
  294.  *   usage: input (argc, argv, &pct, &wSize, &hSize,
  295.  *     &fltrSize, &fltrVar, &threshG,
  296.  *     &gFlag, &qfFlag, &gfFlag)
  297.  */
  298.  
  299.  
  300. #define USAGE_EXIT(VALUE) {usage (VALUE); return (-1);}
  301.  
  302. long
  303. input (argc, argv, pct, wSize, hSize, fltrSize, fltrVarNum, threshG,
  304.        qfFlag, gfFlag)
  305.      int argc;
  306.      char *argv[];
  307.      double *pct;               /* subsample percentage */
  308.      long *wSize, *hSize;       /* subsampled image size: width, height */
  309.      long *fltrSize;            /* filter sidelength */
  310.      double *fltrVarNum;        /* num. of var.s at max width of fltr */
  311.      double *threshG;           /* gray-scale threshold scale factor */
  312.      short *qfFlag;             /* quick filter if =1, or default if 0 */
  313.      short *gfFlag;             /* Gaussian filter if =1, or default if 0 */
  314. {
  315.   long n;
  316.   double atof ();
  317.  
  318.   if (argc < 3)
  319.     USAGE_EXIT (-1);
  320.  
  321.   *pct = DFLTPCT;
  322.   *wSize = *hSize = -1;
  323.   *fltrSize = 1;
  324.   *threshG = -1.0;
  325.   *fltrVarNum = FLTR_VAR_NUM;
  326.   *qfFlag = 0;
  327.   *gfFlag = 0;
  328.  
  329.   for (n = 3; n < argc; n++) {
  330.     if (strcmp (argv[n], "-p") == 0) {
  331.       if (++n == argc || argv[n][0] == '-')
  332.         USAGE_EXIT (-1);
  333.       *pct = atof (argv[n]);
  334.     }
  335.     else if (strcmp (argv[n], "-w") == 0) {
  336.       if (++n == argc || argv[n][0] == '-')
  337.         USAGE_EXIT (-1);
  338.       *wSize = atol (argv[n]);
  339.     }
  340.     else if (strcmp (argv[n], "-h") == 0) {
  341.       if (++n == argc || argv[n][0] == '-')
  342.         USAGE_EXIT (-1);
  343.       *hSize = atol (argv[n]);
  344.     }
  345.     else if (strcmp (argv[n], "-f") == 0) {
  346.       if (++n == argc || argv[n][0] == '-')
  347.         USAGE_EXIT (-1);
  348.       *fltrSize = atol (argv[n]);
  349.     }
  350.     else if (strcmp (argv[n], "-tg") == 0) {
  351.       if (++n == argc || argv[n][0] == '-')
  352.         USAGE_EXIT (-1);
  353.       *threshG = atof (argv[n]);
  354.     }
  355.     else if (strcmp (argv[n], "-v") == 0) {
  356.       if (++n == argc || argv[n][0] == '-')
  357.         USAGE_EXIT (-1);
  358.       *fltrVarNum = atof (argv[n]);
  359.     }
  360.     else if (strcmp (argv[n], "-qf") == 0)
  361.       *qfFlag = 1;
  362.     else if (strcmp (argv[n], "-gf") == 0)
  363.       *gfFlag = 1;
  364.     else if (strcmp (argv[n], "-L") == 0) {
  365.       print_sos_lic ();
  366.       exit (0);
  367.     }
  368.     else
  369.       USAGE_EXIT (-1);
  370.   }
  371.  
  372.   return (0);
  373. }
  374.